#!/bin/csh
#Common stuff for gmt jobs of tisc:

#for old job versions to be compatible with this script:
set name = $prj

#Setting up variables:
set ps  	= $prj.ps
set tmp 	= tmp.tisc.gmt.$prj

if ($prj == "") then 
	echo NEED TO KNOW THE PROJECT NAME!
	exit
endif

if ($?width == "0") set width = 8
if ($?separation == "0") set separation = .8
if ($?intens == "0") set intens	= .2
if ($?height_CS == "0") set height_CS = 3.5

if (-r $prj.hrz) then
	awk '{\
		if (NR==4) {numhrzs=NF-2; ymax=$2; sigx=sigy=0; zmin=0; zmax=0}\
		if (NR>4) {\
			if (zmax<$NF && $NF<10000) zmax=$NF;\
			if (zmin>$3 && $3>-20000)  zmin=$3;\
			if (sigx==1 && $2!=ymin) sigx=0;\
			if (sigx==1) {xmax=$1; countx++;}\
			if ($2!=ymax && sigy==0) {countx = 1; county = 2; xmin=$1; ymin=$2; sigx=sigy=1;}\
			if ($2!=ymin) {ymin=$2; county++;}\
		}\
	}\
	END {print numhrzs, countx, county, xmin, xmax, ymin, ymax, zmin/1e3, zmax/1e3}' $prj.hrz > $tmp.prop.tmp
else
	echo ERROR: File $prj.hrz does not exist\!\!
	exit
endif

set props = `cat $tmp.prop.tmp` 
set numpfls = $props[1]
set Nx = $props[2]
set Ny = $props[3]
set x0 = $props[4]
set xf = $props[5]
set y0 = $props[6]
set yf = $props[7]
set Region = $x0/$xf/$y0/$yf

set x0m = `echo  $x0 \* 1000 | bc -l`
set xfm = `echo  $xf \* 1000 | bc -l`
set y0m = `echo  $y0 \* 1000 | bc -l`
set yfm = `echo  $yf \* 1000 | bc -l`

set Lx	= `echo  $xf - $x0  | bc -l`
set Ly	= `echo  $yf - $y0  | bc -l`
set dx	= `echo $Lx \/ \($Nx-1\) | bc -l`
set dy	= `echo $Ly \/ \($Ny-1\) | bc -l`
set dxy = `awk -v dx=$dx -v dy=$dy 'BEGIN{print sqrt(dx*dx+dy*dy)}'`
if (-r $prj.PRM) then
	set dt   	= `awk '{if ($1=="dt")        {print $2; exit;}}' $prj.PRM`
	set Timeini 	= `awk '{if ($1=="Timeini")   {print $2; exit;}}' $prj.PRM`
	set Timefinal 	= `awk '{if ($1=="Timefinal") {print $2; exit;}}' $prj.PRM`
	set zmin 	= `awk '{if ($1=="zmin")      {print $2/1e3; exit;}}' $prj.PRM`
	set zmax 	= `awk '{if ($1=="zmax")      {print $2/1e3; exit;}}' $prj.PRM`
endif
if ($zmin == "") set zmin = $props[8]
if ($zmax == "") set zmax = $props[9]
set zminm = `echo  $zmin \* 1000 | bc -l`
set zmaxm = `echo  $zmax \* 1000 | bc -l`
set ZRegion	= $zmin/$zmax
set Timenow = `awk -v dt=$dt '(NR==1) {txt=substr($(NF-1),4,8); if (dt>.1) printf("%.1f", txt); else printf("%.2f", txt); exit;}' $prj.hrz`
set xtime = `echo $x0 + $Lx \* .001  | bc -l`
set ytime = `echo $yf - $Ly \* .001  | bc -l`
set xtitle = `echo $xf - $Lx \* .007  | bc -l`
set ytitle = `echo $y0 + $Ly \* .015 | bc -l`

set half_x = `echo $width + \( $separation \) \/ 2 | bc -l`
set half_separation = `echo $separation \/ 2 | bc -l`
set half_width = `echo $width \/ 2 | bc -l`
set total_width = `echo $width \* 2 + $separation | bc -l`
set height = `echo $width \* $Ly \/ $Lx  | bc -l`
set size = $width/$height
set vert_shift = `echo -$separation - $height  | bc -l`
set horz_shift = `echo  $separation + $width   | bc -l`
set xt = `echo $xf - 12  | bc -l`
set yt = `echo $y0 + 12  | bc -l`
set vert_shift_CS = `echo -$height_CS \- $separation | bc -l`
set width_CS = $width
set horz_shift_CS = 0
if (-r $prj.pfl) then 
	set dmax = `tail $prj.pfl | awk '{dmax=$3;} END{print dmax}'`
	set dens = `awk '(NR==2) {print $0; exit;}' $prj.pfl`
	set age  = `awk '(NR==3) {print $0; exit;}' $prj.pfl`
	set sedim_age_limit = `echo \( $Timeini + $Timefinal \) / 2 | bc -l` 
endif


#echo
#echo -n Region = $Region   Nx=$Nx  Ny=$Ny

if ($2 != "" && $5 != "") then
	set Region_plot = $2/$3/$4/$5 
	echo Region_plot = $Region_plot  ZRegion=$ZRegion
	set Lx_plot	= `echo  $3 - $2  | bc -l`
	set Ly_plot	= `echo  $5 - $4  | bc -l`
else 
	set Region_plot = $Region
	set Lx_plot	= $Lx
	set Ly_plot	= $Ly
endif

#size of rivers depending on domain area:
set area = `echo  $Lx_plot \* $Ly_plot | bc -l`
set meanL = `echo  \( $Lx_plot \+ $Ly_plot \) \/ 2 | bc -l`
set Lxy = `awk -v Lx=$Lx_plot -v Ly=$Ly_plot 'BEGIN{print sqrt(Lx*Lx+Ly*Ly)}'`
set disch1 = `echo  $area \/ 90000 \* 1    | bc -l`
set disch2 = `echo  $area \/ 90000 \* 10   | bc -l`
set disch3 = `echo  $area \/ 90000 \* 100  | bc -l`
set sedld1 = `echo  $area \/ 90000 \* .20  | bc -l`
set sedld2 = `echo  $area \/ 90000 \* 2.0  | bc -l`
set sedld3 = `echo  $area \/ 90000 \* 20   | bc -l`
set sedloadmax = `awk -v sedld3=$sedld3 'BEGIN{print sedld3*5}' `  #`awk '{print $4/$3}' $tmp.chosen_basin.bas.tmp | sort -n | tail -n 1`

set vert_scale	= 1.5
set view_ang	= 210/30

set color_heavy_basement = 140/100/0
set color_basement = 165/120/0
set color_heavy_sediment = 80/130/255
set color_old_sediment = 255/180/0
set color_sediment = 255/255/20
set color_basement = 165/120/0
set color_thrust   = 170/115/10
set color_sea      = 50/180/255
set col_lake 	= 0/40/255
set col_river	= 0/0/255
set color_ice 	= 220/220/255
set col_lin	= 4/255/0/120
set col_masstr	= 0/0/0

set ticklabels 	= `awk -v meanL=$meanL 'BEGIN{if (meanL>=990) print "a500f50"; else if (meanL>=390) print "a200f50"; else if (meanL>=90) print "a50f10"; else print "a10f1";}'`
set zticklabels	= `awk -v zmin=$zmin -v zmax=$zmax 'BEGIN{if (zmax-zmin>=30) print "a10f5"; else if (zmax-zmin>=10) print "a4f1"; else if (zmax-zmin>=3) print "a1f.5"; else print "a10f1";}'`

if (-r $prj.xyw) then
	set sea_level = `awk '{if (NR==1) printf ("%.0f",$(NF-1)); exit;}' $prj.xyw`
	set ndrain = `awk '{if ($5 != "-") n++}END{print n}' $prj.xyw`
else
	set sea_level = 0
	set ndrain = 0
endif


#CORNERS IN PRFL FILE:
if (-r $prj.PRFL) then
set x_corner = ( 1 2 3 4 5 6 7 8 9 10 11 12 13 )
set x_corner[1] = `awk '{a=$1; b=$2; if (NR>1) {l += sqrt((a-aa)*(a-aa)+(b-bb)*(b-bb)); if (NR==2) print l/1e3}; aa=a; bb=b}' $prj.PRFL`
set x_corner[2] = `awk '{a=$1; b=$2; if (NR>1) {l += sqrt((a-aa)*(a-aa)+(b-bb)*(b-bb)); if (NR==3) print l/1e3}; aa=a; bb=b}' $prj.PRFL`
set x_corner[3] = `awk '{a=$1; b=$2; if (NR>1) {l += sqrt((a-aa)*(a-aa)+(b-bb)*(b-bb)); if (NR==4) print l/1e3}; aa=a; bb=b}' $prj.PRFL`
set x_corner[4] = `awk '{a=$1; b=$2; if (NR>1) {l += sqrt((a-aa)*(a-aa)+(b-bb)*(b-bb)); if (NR==5) print l/1e3}; aa=a; bb=b}' $prj.PRFL`
set x_corner[5] = `awk '{a=$1; b=$2; if (NR>1) {l += sqrt((a-aa)*(a-aa)+(b-bb)*(b-bb)); if (NR==6) print l/1e3}; aa=a; bb=b}' $prj.PRFL`
set x_corner[6] = `awk '{a=$1; b=$2; if (NR>1) {l += sqrt((a-aa)*(a-aa)+(b-bb)*(b-bb)); if (NR==7) print l/1e3}; aa=a; bb=b}' $prj.PRFL`
set x_corner[7] = `awk '{a=$1; b=$2; if (NR>1) {l += sqrt((a-aa)*(a-aa)+(b-bb)*(b-bb)); if (NR==8) print l/1e3}; aa=a; bb=b}' $prj.PRFL`
set x_corner[8] = `awk '{a=$1; b=$2; if (NR>1) {l += sqrt((a-aa)*(a-aa)+(b-bb)*(b-bb)); if (NR==9) print l/1e3}; aa=a; bb=b}' $prj.PRFL`
set x_corner[9] = `awk '{a=$1; b=$2; if (NR>1) {l += sqrt((a-aa)*(a-aa)+(b-bb)*(b-bb)); if (NR==10) print l/1e3}; aa=a; bb=b}' $prj.PRFL`
set x_corner[10] = `awk '{a=$1; b=$2; if (NR>1) {l += sqrt((a-aa)*(a-aa)+(b-bb)*(b-bb)); if (NR==11) print l/1e3}; aa=a; bb=b}' $prj.PRFL`
set x_corner[11] = `awk '{a=$1; b=$2; if (NR>1) {l += sqrt((a-aa)*(a-aa)+(b-bb)*(b-bb)); if (NR==12) print l/1e3}; aa=a; bb=b}' $prj.PRFL`
set x_corner[12] = `awk '{a=$1; b=$2; if (NR>1) {l += sqrt((a-aa)*(a-aa)+(b-bb)*(b-bb)); if (NR==13) print l/1e3}; aa=a; bb=b}' $prj.PRFL`
set x_corner[13] = `awk '{a=$1; b=$2; if (NR>1) {l += sqrt((a-aa)*(a-aa)+(b-bb)*(b-bb)); if (NR==14) print l/1e3}; aa=a; bb=b}' $prj.PRFL`
#echo corners: $x_corner[1] $x_corner[2] $x_corner[11]
endif

#PALLETTES:
generate_topo_cpt.job $sea_level > $tmp.topo_palet.tmp
#makecpt -Cglobe  > $tmp.topo_palet.tmp


cat - <<END>  $tmp.ice_cpt.tmp
	0	210	210	235		200	210	210	235
	200	230	210	255		1000	230	210	255
B	0	0	0
F	255	255	255
END

cat - <<END>  $tmp.lakes_cpt.tmp
	0	0	0	0	.49	0	0	0
	.49	0	120	255	1	0	120	255
END
cat - <<END>  $tmp.lakes_endorh_cpt.tmp
	0	0	0	0	.49	0	0	0
	.49	30	150	255	1	30	150	255
END

set min_rain = .5
cat - <<END>  $tmp.rain_cpt.tmp 
#0		255	150	0	$min_rain 255	255	0
$min_rain 	21      234     255     1     21      234     255
1        	64      200     255     1.5       64      200     255
1.5       	106     149     255     2     106     149     255
2       	149     106     255     3       149     106     255
3       	212     43      255     5       212     43      255
B       0       0       0
F       230     0       255
N       128     128     190
END

cat - <<END>  $tmp.subs_cpt.tmp
	-5000	255	0	0	-500	255	80	80
	-500	255	80	80	-50	255	170	170
	-50	255	170	170	0	255	255	255
	0	255	255	255	200	170	170	255
	200	170	170	255	2000	80	80	255
	2000	80	80	255	20000	0	0	255
B	100 0 0
F	0 0 60
END

cat - <<END> $tmp.deflect_cpt.tmp
-2	0	255	0	0	100	255	200
 0	150	255	255	2	64	159	255
 2	64	159	255	4	64	127	255
 4	64	127	255	6	192	0	0
B	159	180	64
F	192	0	0
END

cat - <<END>  $tmp.subsrate_cpt.tmp
	-1000	200	0	0	-300	255	60	180
	-300	255	60	180	-100	255	160	60
	-100	255	160	60	0	255	255	255
	0	255	255	255	100	180	255	60
	100	180	255	60	200	60	255	160
	200	60	255	160	400	0	200	0
B	190 0 0
F	0 190 0
END

cat - <<END>  $tmp.eet_cpt.tmp
	0	255	255	255	5	255	255	0
	5	255	255	0	10	255	0	0
	10	255	0	0	30	0	0	255
	30	0	0	255	80	0	0	0
END

cat - <<END>  $tmp.sed_cpt.tmp
	0	255	255	255	.2	210	255	170
	.2	210	255	170	1	255	255	0
	1	255	255	0	2	255	255	0
	2	255	190	0	4	255	190	0
	4	255	120	0	6	255	120	0
	6	255	0	0	8	255	0	0
	8	155	0	0	10	155	0	0
B	255 255 255
F	30   30   30
END

cat - <<END>  $tmp.erosrate_cpt.tmp
	-10	255	0	0	-1	150	190	0
	-1	150	190	0	-.1	255	220	30
	-.1	255	220	30	0	200	200	200
	0	200	200	200	.02	0	200	255
	.02	0	200	255	.1	90	160	255
	.1	128	127	255	1	128	127	255
	1	212	43	255	10	212	43	255
END
#makecpt -Ccool -T.1/10/1 -Qo -D >> $tmp.erosrate_cpt.tmp

cat - <<END>  $tmp.eros_cumul_cpt.tmp
	-8	150	0	200	-3	150	0	200
	-3	70	0	255	-2	70	0	255
	-2	0	60	255	-1	0	60	255
	-1	30	140	255	0	255	255	255
END
cat $tmp.sed_cpt.tmp >>  $tmp.eros_cumul_cpt.tmp

#makecpt -T0/$sedloadmax/1 > $tmp.sedload_cpt.tmp
cat  <<END> $tmp.sedload_cpt.tmp
	0       30   180     255     		.1          200   200    0
	.1    200   200     0     		.2           255   100    0
	.2    255   100     0     		0.5         255    00    0
	0.5   255    00     0     		$sedloadmax       50     0    0
F	0	0	0
B	255	255	255
END

cat - <<END>  $tmp.geol_palet.tmp
	500	40	80	255	1500	40	80	255
	1500	255	255	0	2045	255	255	0
	2045	255	190	0	2399	255	190	0
	2399	210	120	80	2775	210	120	80
	2775	165	120	0	2825	165	120	0
	2825	155	120	0	3500	155	120	0
B	30 30 30
F	255 0 0
END

if (0) then 
  #BLACK AND WHITE B/W pallette
  set color_basement = 100/100/100
  set color_thrust   = 120/120/120
  set color_sediment = 200/200/200
  set color_sea      = 255/255/255
  set col_lin	  = 6/0
  set col_river   = 0/0/0
  set col_masstr  = 0/0/0
  cat - <<END>  $tmp.geol_palet.tmp
  	  500	  200	  200	  255	  1500    200	  200	  255
  	  1500    200	  200	  200	  2500    200	  200	  200
  	  2500    120	  120	  120	  2825    120	  120	  120
  	  2825    100	  100	  100	  3500    100	  100	  100
  B	  30 30 30
  F	  255 0 0
END
endif


#GRIDS:
#Topography
if (-r $prj.xyw) then 
awk -v sl=$sea_level '{topo=$6; if ($5 =="S" && topo>=sl) topo=sl-1.11; print($1,$2,topo)}' $prj.xyw |\
	xyz2grd	-G$tmp.top.grd.tmp -I$dx/$dy -R$Region -H2
else
awk '{print($1,$2,$NF)}' $prj.hrz |\
	xyz2grd	-G$tmp.top.grd.tmp -I$dx/$dy -R$Region -H2
endif

#Ilumination:
grdgradient $tmp.top.grd.tmp -A190 -G$tmp.intensity.top.grd.tmp  #-Nt1
#scales with model length. 00002 MUL is normal for orogen+basin scale (1000km)
grdmath $tmp.intensity.top.grd.tmp $meanL MUL .00001 MUL = $tmp.intensity.top_all.grd.tmp
#iluminate only above sea
grdmath $tmp.top.grd.tmp 0 GT = $tmp.contin.grd.tmp
grdmath $tmp.intensity.top_all.grd.tmp $tmp.contin.grd.tmp MUL = $tmp.intensity.top.grd.tmp
#grdinfo $tmp.intensity.top.grd.tmp

#Lakes:
if (-r $prj.lak) then 
	awk -v sl=$sea_level '(substr($0,1,1)!="#"){if (substr($0,1,1)==">") disch=$11; else{if (disch>dl && ($3=="L" || $3=="E")) print $1,$2,1; else print $1,$2,0}}' dl=$disch1 $prj.lak | \
		xyz2grd -G$tmp.lakes.grd.tmp -I$dx/$dy -R$Region -N0
	#grdcontour $tmp.lakes.grd.tmp -JX1 -R -C$tmp.lakes_cpt.tmp -M -D$tmp.lakes.xy.tmp -O -K >> $ps
	awk -v sl=$sea_level '(substr($0,1,1)!="#"){if (substr($0,1,1)==">") {disch=$11; outl=$9;} else{if (disch>dl && ($3=="L" || $3=="E") && outl==0) print $1,$2,1; else print $1,$2,0}}' dl=$disch1 $prj.lak | \
		xyz2grd -G$tmp.lakes_endorh.grd.tmp -I$dx/$dy -R$Region -N0
	#grdcontour $tmp.lakes_endorh.grd.tmp -JX1 -R -C$tmp.lakes_cpt.tmp -M -D$tmp.lakes_endorh.xy.tmp -O -K >> $ps
endif

#Sediment thickness:
if (-r $prj.hrz) then 
	awk '{if (NR==2) for (i=3; i<=NF; i++) {dens[i]=$i;}\
	    if (NR>2) { sedthick=0; \
		for (i=4; i<=NF; i++) {if (dens[i]<2500) sedthick+=($i-$(i-1)); }\
		print($1,$2,sedthick/1e3); \
	} }' $prj.hrz | xyz2grd -G$tmp.sed.grd.tmp -I$dx/$dy -R -H0
endif

#eros./sed.:
if (-r $prj.st) then
	awk '{print($1,$2, $4/1e3);}' $prj.st | xyz2grd -G$tmp.eros_cumul.grd.tmp -I$dx/$dy -R -H2
	awk '{print($1,$2, $5/1e3);}' $prj.st | xyz2grd -G$tmp.erosrate.grd.tmp   -I$dx/$dy -R -H2
endif

#Lithology/Geology:
if (-r $prj.hrz) then
	awk '{if (NR==2) for (i=3; i<=NF; i++) dens[i]=$i; if (NR>2) {for (i=NF; i>3; i--) if ($i-$(i-1)>1) break; if ($NF>=0) print $1,$2,dens[i]; else print $1,$2,1000}}' $prj.hrz > $tmp.geol.xyd.tmp
	xyz2grd	$tmp.geol.xyd.tmp -G$tmp.geol.grd.tmp -I$dx/$dy -R$Region
endif
